#ifndef _BERNSTEIN_H_
#define _BERNSTEIN_H_
__device__ void GetBernstein3(float Bern[], float u);

//__device__ void GetBernstein2(float Bern[], float u);


__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re);

//__device__ void EvalBezierBlossomNormal2(bool isOne, float u, float v, int edgei, int ui, int vi, int nx, int ny, int nz, float* func, float* re, float* norm);
//
//__device__ void EvalBezierBlossomPolarTex(float u, float v, float w, int nx, int ny, int nz, int index, float* re);
//
//
//__device__ void EvalBezierBlossomNormal(float u, float v, float w, int nx, int ny, int nz, float* func, float* re);
//
//__device__ void EvalBezierBlossomNormal2deg(float u, float v, float w, float* func, float* re, float* norm);

//__device__ void GetBernsteinfw(float Bern[], int deg, float u);
//__device__ void EvalBezierBlossomPolarTexfw(float u, float v, float w, int nx, int ny, int nz, int index, float* re);
//__device__ void EvalBezierBlossomTex(float u, float v, float w, int nx, int ny, int nz, texture<float, 1, cudaReadModeElementType> tx, int index, float* re);

__device__ void GetBernstein3(float Bern[], float u)
{
	float cu = 1.0f-u;
	Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u;
}

__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}
__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	int o = 0;
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += func[o++] * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

//__device__ void GetBernstein2(float Bern[], float u)
//{
//	float cu = 1.0f-u;
//	Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u;
//}


__device__ void EvalBezierBlossom(float u, float v, float w, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein3(nn, w);
	for (int i=0; i<=3; i++)
	for (int j=0; j<=3; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=3; k++)
			f1[i][j] += *func++ * nn[k];
	}

	GetBernstein3(nn, v);
	for (int i=0; i<=3; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=3; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein3(nn, u);
	(*re) = 0.0f;
	for (int i=0; i<= 3; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossom3(float u, float v, float w, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein3(nn, w);
	f1[0][0] = func[ 0]*nn[0]+func[ 1]*nn[1]+func[ 2]*nn[2]+func[ 3]*nn[3];
	f1[0][1] = func[ 4]*nn[0]+func[ 5]*nn[1]+func[ 6]*nn[2]+func[ 7]*nn[3];
	f1[0][2] = func[ 8]*nn[0]+func[ 9]*nn[1]+func[10]*nn[2]+func[11]*nn[3];
	f1[0][3] = func[12]*nn[0]+func[13]*nn[1]+func[14]*nn[2]+func[15]*nn[3];
	f1[1][0] = func[16]*nn[0]+func[17]*nn[1]+func[18]*nn[2]+func[19]*nn[3];
	f1[1][1] = func[20]*nn[0]+func[21]*nn[1]+func[22]*nn[2]+func[23]*nn[3];
	f1[1][2] = func[24]*nn[0]+func[25]*nn[1]+func[26]*nn[2]+func[27]*nn[3];
	f1[1][3] = func[28]*nn[0]+func[29]*nn[1]+func[30]*nn[2]+func[31]*nn[3];
	f1[2][0] = func[32]*nn[0]+func[33]*nn[1]+func[34]*nn[2]+func[35]*nn[3];
	f1[2][1] = func[36]*nn[0]+func[37]*nn[1]+func[38]*nn[2]+func[39]*nn[3];
	f1[2][2] = func[40]*nn[0]+func[41]*nn[1]+func[42]*nn[2]+func[43]*nn[3];
	f1[2][3] = func[44]*nn[0]+func[45]*nn[1]+func[46]*nn[2]+func[47]*nn[3];
	f1[3][0] = func[48]*nn[0]+func[49]*nn[1]+func[50]*nn[2]+func[51]*nn[3];
	f1[3][1] = func[52]*nn[0]+func[53]*nn[1]+func[54]*nn[2]+func[55]*nn[3];
	f1[3][2] = func[56]*nn[0]+func[57]*nn[1]+func[58]*nn[2]+func[59]*nn[3];
	f1[3][3] = func[60]*nn[0]+func[61]*nn[1]+func[62]*nn[2]+func[63]*nn[3];
	//f1[0][0] = 0;
	//f1[0][1] = 0;
	//f1[0][2] = 0;
	//f1[0][3] = 0;
	//f1[1][0] = 0;
	//f1[1][1] = 0;
	//f1[1][2] = 0;
	//f1[1][3] = 0;
	//f1[2][0] = 0;
	//f1[2][1] = 0;
	//f1[2][2] = 0;
	//f1[2][3] = 0;
	//f1[3][0] = 0;
	//f1[3][1] = 0;
	//f1[3][2] = 0;
	//f1[3][3] = 0;

	GetBernstein3(nn, v);
	f2[0] = f1[0][0]*nn[0]+f1[0][1]*nn[1]+f1[0][2]*nn[2]+f1[0][3]*nn[3];
	f2[1] = f1[1][0]*nn[0]+f1[1][1]*nn[1]+f1[1][2]*nn[2]+f1[1][3]*nn[3];
	f2[2] = f1[2][0]*nn[0]+f1[2][1]*nn[1]+f1[2][2]*nn[2]+f1[2][3]*nn[3];
	f2[3] = f1[3][0]*nn[0]+f1[3][1]*nn[1]+f1[3][2]*nn[2]+f1[3][3]*nn[3];

	GetBernstein3(nn, u);
	(*re) = f2[0]*nn[0]+f2[1]*nn[1]+f2[2]*nn[2]+f2[3]*nn[3];

}

//__device__ void EvalBezierBlossomNormal2(bool isOne, float u, float v, int edgei, int ui, int vi, int nx, int ny, int nz, float* func, float* re, float* norm)
//{
//	float nuv[4], tempv[4];
//	float norm1,norm2;
//
//	float nn[4];
//	float temp[4];
//	float deg[3] = {nx, ny, nz};
//	int ind[3] = {0,0,0};
//	float funcInd=0.0f;
//
//	if(isOne)
//		ind[edgei] = deg[edgei];
//	else
//		ind[edgei] = 0;
//
//
//	GetBernstein(nn, deg[vi], v);
//	GetBernstein(nuv, deg[vi], v+DERIVATIVE_DELTA);
//	for(int i=0; i<=deg[ui]; i++)
//	{
//		temp[i] = 0.0f;
//		tempv[i] = 0.0f;
//		ind[ui] = i;
//		for(int j=0; j<=deg[vi]; j++)
//		{
//			ind[vi] = j;
//			funcInd = *(func+(ind[0]*(ny+1)+ind[1])*(nz+1)+ind[2]);
//			temp[i] += funcInd*nn[j];
//			tempv[i] += funcInd*nuv[j];
//		}
//	}
//
//
//
//	GetBernstein(nn, deg[ui], u);
//	GetBernstein(nuv, deg[ui], u+DERIVATIVE_DELTA);
//	(*re) = 0.0f;
//	norm1 = 0.0f;
//	norm2 = 0.0f;
//	for(int i=0; i<=deg[ui]; i++)
//	{
//		(*re) += temp[i]*nn[i];
//		norm1 += temp[i]*nuv[i];
//		norm2 += tempv[i]*nn[i];
//	}
//	(*norm) = (norm1-(*re))*DERIVATIVE_DELTA_INV;
//	(*(norm+1)) = (norm2-(*re))*DERIVATIVE_DELTA_INV;
//}
//
//
__device__ void EvalBezierBlossomNormal(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}
	//func += degShift[0];
	func += (ny+1)*(nz+1);
	GetBernstein(nn, ny, v);
	for (int i=0; i<nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx-1, u);
	re[0] = 0.0f;
	for (int i=0; i<nx; i++)
	{
		re[0] += f2[i]*nn[i];
	}

	//Y direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<=nz; k++)
				f1[i][j] += *func++ * nn[k];
		}
		//func += degShift[1];
		func += (nz+1);
	}
	GetBernstein(nn, ny-1, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[1] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[1] += f2[i]*nn[i];
	}

	//Z direction:
	GetBernstein(nn, nz-1, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<=ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<nz; k++)
				f1[i][j] += *func++ * nn[k];
			func++;
		}
	}
	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[2] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[2] += f2[i]*nn[i];
	}
}
//
//__device__ void EvalBezierBlossomNormal2deg(float u, float v, float w, float* func, float* re, float* norm)
//{
//	float nn[3], deri[3];
//	GetBernstein2(nn, w);
//	GetBernstein2(deri, w+DERIVATIVE_DELTA);
//	float coeff=0.0f;
//	
//	/*********************************************************************/
//	float f1[9], z1[9];
//	for(int i=0; i<9; i++)
//	{		
//		coeff = *func++;
//		f1[i] = coeff * nn[0];
//		z1[i] = coeff * deri[0];
//		coeff = *func++;
//		f1[i] += coeff * nn[1];
//		z1[i] += coeff * deri[1];
//		coeff = *func++;
//		f1[i] += coeff * nn[2];
//		z1[i] += coeff * deri[2];
//	}
//	
//	float f2[3], z2[3], y2[3];
//	GetBernstein2(nn, v);
//	GetBernstein2(deri, v+DERIVATIVE_DELTA);
//	f2[0] = f1[0]*nn[0]+f1[1]*nn[1]+f1[2]*nn[2];
//	f2[1] = f1[3]*nn[0]+f1[4]*nn[1]+f1[5]*nn[2];
//	f2[2] = f1[6]*nn[0]+f1[7]*nn[1]+f1[8]*nn[2];	
//		
//	z2[0] = z1[0]*nn[0]+z1[1]*nn[1]+z1[2]*nn[2];
//	z2[1] = z1[3]*nn[0]+z1[4]*nn[1]+z1[5]*nn[2];
//	z2[2] = z1[6]*nn[0]+z1[7]*nn[1]+z1[8]*nn[2];	
//		
//	y2[0] = f1[0]*deri[0]+f1[1]*deri[1]+f1[2]*deri[2];
//	y2[1] = f1[3]*deri[0]+f1[4]*deri[1]+f1[5]*deri[2];
//	y2[2] = f1[6]*deri[0]+f1[7]*deri[1]+f1[8]*deri[2];	
//	/*********************************************************************/
//	
//
//	////These commened lines can replace the codes surround by "****" line.
//	//float f1[3][3], z1[3][3];
//	//for (int i=0; i<=2; i++)
//	//for (int j=0; j<=2; j++)
//	//{
//	//	f1[i][j] = 0.0f;
//	//	z1[i][j] = 0.0f;
//	//	for (int k=0; k<=2; k++)
//	//	{
//	//		coeff = *func++;
//	//		f1[i][j] += coeff * nn[k];
//	//		z1[i][j] += coeff * deri[k];
//	//	}
//	//}
//
//	//
//	//float f2[3], z2[3], y2[3];
//	//GetBernstein2(nn, v);
//	//GetBernstein2(deri, v+DERIVATIVE_DELTA);
//	//for (int i=0; i<=2; i++)
//	//{
//	//	f2[i] = 0.0f;
//	//	y2[i] = 0.0f;
//	//	z2[i] = 0.0f;
//	//	for (int j=0; j<=2; j++)
//	//	{
//	//		f2[i] += f1[i][j]*nn[j];
//	//		y2[i] += f1[i][j]*deri[j];
//	//		z2[i] += z1[i][j]*nn[j];
//	//	}
//	//}
//
//
//	GetBernstein2(nn, u);
//	float fValue = f2[0]*nn[0]+f2[1]*nn[1]+f2[2]*nn[2];
//	norm[1] = y2[0]*nn[0]+y2[1]*nn[1]+y2[2]*nn[2];
//	norm[2] = z2[0]*nn[0]+z2[1]*nn[1]+z2[2]*nn[2];
//	
//	GetBernstein2(deri, u+DERIVATIVE_DELTA);
//	norm[0] = f2[0]*deri[0]+f2[1]*deri[1]+f2[2]*deri[2];
//	
//	norm[0] = (norm[0]-fValue)*DERIVATIVE_DELTA_INV;
//	norm[1] = (norm[1]-fValue)*DERIVATIVE_DELTA_INV;
//	norm[2] = (norm[2]-fValue)*DERIVATIVE_DELTA_INV;
//	(*re) = fValue;
//}
//
//__device__ void GetIntBernstein(int Bern[], int deg, int part, int u)
//{
//	//Multiplication of int will be overflow!
//	int cu = part-u;
//
//	switch (deg)
//	{
//	case 2:
//		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
//	case 3:
//		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
//	case 1:
//		Bern[0]=cu; Bern[1]=u; break;
//	}
//}
//
//__device__ void GetCoeffPartBern(int coeffO, int ind, int* partB1, int* partB2, int part, int deg)
//{
//	int part2 = ind%(part*part);
//	int part1 = part2/part;
//	part2 = part2%part;
//
//	GetIntBernstein(partB1, deg, part, part1);
//	GetIntBernstein(partB2, deg, part, part2);
//}
//


#endif
